rm(list = ls())
library(ggplot2)
library(matlab)
library(patchwork)
# parameters setting
psi <- 0.05
theta <- 0.5
beta1 <- 20
beta2 <- 0.1
mu <- 0.01

# exo variables
m0 <- 100
beta0 <- 500
ytn <- 2000
p0star <- 0
i0star <- 3

ode_solve <- function(psi, theta, beta1,beta2, mu, m0, beta0,ytn,p0star,i0star){
  # coef matrix
  A <- matrix(c(-mu*(beta1 + beta2/theta),mu*beta1,1/theta,0),2,2,byrow = T)
  B <- matrix(c(mu,mu*beta2/theta,-mu*(psi*beta2/theta+1),0,mu*beta1,0,-1/theta,psi/theta,-1,0),nrow = 2,byrow = T)
  z <- matrix(c(beta0,m0,ytn,i0star,p0star),ncol = 1)

  # steady state
  st <- -solve(A) %*% B %*% z
  # eig value
  ev <- eigen(A)$values
  return(list(A=A, B=B, z=z, st=st, ev=ev))
}

# 稳态
rlt <- ode_solve(psi, theta, beta1,beta2, mu, m0, beta0,ytn,p0star,i0star)
rlt$st

# 特征值
rlt$ev

# 脉冲响应分析
# 新的稳态
m0 <- 101
rlt_new <- ode_solve(psi, theta, beta1,beta2, mu, m0, beta0,ytn,p0star,i0star)
st_new <- rlt_new$st

# 冲击下，第一期的值
jmp <- (rlt_new$A + eye(2)) %*% rlt$st + rlt_new$B %*% rlt_new$z
# 汇率的瞬间跳跃
jmp[2] <- -(-rlt_new$ev[1]*theta*st_new[2]+i0star*theta-psi*ytn+m0-jmp[1])/(rlt_new$ev[1]*theta)

# 鞍点路径的运动
picdata <- data.frame(p = c(rlt$st[1],jmp[1]),s = c(rlt$st[2],jmp[2]))
for (i in 3:15) {
  picdata[i,] <- 1/(1-rlt_new$ev[1])*picdata[i-1,]-rlt_new$ev[1]/(1-rlt_new$ev[1])*st_new
}
picdata$i <- -(m0-picdata$p-psi*ytn)/theta
picdata$i[1] <- -(100-picdata$p[1]-psi*ytn)/theta # 修正初期值,即用原来的m0

picdata$ytd <- beta0 + beta1*(picdata$s-picdata$p+p0star)+beta2/theta*(m0-picdata$p-psi*ytn)
picdata$ytd[1] <- beta0 + beta1*(picdata$s[1]-picdata$p[1]+p0star)+beta2/theta*(100-picdata$p[1]-psi*ytn)
picdata$time <- 1:nrow(picdata)

pexc <- ggplot(picdata, aes(x = time, y = s)) + geom_line() + theme_bw()
pp <- ggplot(picdata, aes(x = time, y = p)) + geom_line() + theme_bw()
py <- ggplot(picdata, aes(x = time, y = ytd)) + geom_line() + theme_bw()
pi <- ggplot(picdata, aes(x = time, y = i)) + geom_line() + theme_bw()
(pexc / pp)|(py/pi)
# ggsave('../exc.pdf')


# 敏感性分析
mu <- 0.001
rlt <- ode_solve(psi, theta, beta1,beta2, mu, m0=100, beta0,ytn,p0star,i0star)

# 冲击下新的响应路径
rlt_new <- ode_solve(psi, theta, beta1,beta2, mu, m0=101, beta0,ytn,p0star,i0star)
st_new <- rlt_new$st

# 冲击下，第一期的值
jmp <- (rlt_new$A + eye(2)) %*% rlt$st + rlt_new$B %*% rlt_new$z
# 汇率的瞬间跳跃
jmp[2] <- -(-rlt_new$ev[1]*theta*st_new[2]+i0star*theta-psi*ytn+m0-jmp[1])/(rlt_new$ev[1]*theta)

# 鞍点路径的运动
picdata <- data.frame(p = c(rlt$st[1],jmp[1]),s = c(rlt$st[2],jmp[2]))
for (i in 3:15) {
  picdata[i,] <- 1/(1-rlt_new$ev[1])*picdata[i-1,]-rlt_new$ev[1]/(1-rlt_new$ev[1])*st_new
}
picdata$i <- -(m0-picdata$p-psi*ytn)/theta
picdata$i[1] <- -(100-picdata$p[1]-psi*ytn)/theta # 修正初期值,即用原来的m0

picdata$ytd <- beta0 + beta1*(picdata$s-picdata$p+p0star)+beta2/theta*(m0-picdata$p-psi*ytn)
picdata$ytd[1] <- beta0 + beta1*(picdata$s[1]-picdata$p[1]+p0star)+beta2/theta*(100-picdata$p[1]-psi*ytn)
picdata$time <- 1:nrow(picdata)

pexc <- ggplot(picdata, aes(x = time, y = s)) + geom_line() + theme_bw()
pp <- ggplot(picdata, aes(x = time, y = p)) + geom_line() + theme_bw()
py <- ggplot(picdata, aes(x = time, y = ytd)) + geom_line() + theme_bw()
pi <- ggplot(picdata, aes(x = time, y = i)) + geom_line() + theme_bw()
(pexc / pp)|(py/pi)
# `ggsave('../exc_sen.pdf')
